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Abstract 


A cohesive element for shell analysis is presented. The element can be used to simulate 
the initiation and growth of delaminations between stacked, non-coincident layers of 
shell elements. The procedure to construct the element accounts for the thickness offset 
by applying the kinematic relations of shell deformation to transform the stiffness and 
internal force of a zero-thickness cohesive element such that interfacial continuity 
between the layers is enforced. The procedure is demonstrated by simulating the 
response and failure of the Mixed Mode Bending test and a skin-stiffener debond 
specimen. In addition, it is shown that stacks of shell elements can be used to create 
effective models to predict the inplane and delamination failure modes of thick 
components. The results indicate that simple shell models can retain many of the 
necessaiy predictive attributes of much more complex 3D models while providing the 
computational efficiency that is necessaiy for design. 

Keywords: cohesive elements; shell elements; delamination; composite laminates 
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crack length 

area of integration of cohesive element 

strain-displacement matrix 

constitutive matrix 

damage matrix 

damage variable 

length of lever arm in MMB test 

applied force 

internal force of shell cohesive element 

internal force of zero-thickness cohesive element 

mode I, II, and III energy release rates 

critical energy release rate 

critical energy release rates in mode I, II 

total energy release rate 

identity matrix 

J-integral 

penalty stiffness for cohesive law 
stiffness matrix of shell cohesive element 
stiffness matrix of zero-thickness cohesive element 
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element length 
shell transformation matrices 
stress tensor 
pseudo time factor 

pseudo time factor at the maximum load 
thickness of j-th ply 

thickness of the plies above and below an interface 

vector of degrees of freedom at a node for 3D and shell cohesive elements 

vector of all degrees of freedom of 3D and shell cohesive elements 
translational and rotational nodal degrees of freedom 

coefficients of thermal expansion 
displacement jump for damage initiation 
final displacement jump 

mode I, mode II, and mode III displacement jumps 
tangential displacement jump 

difference between the curing temperature and room temperature 
Kronecker delta 

equivalent displacement jump 
viscoelastic regularization factor 

parameter of Benzeggagh-Kenane mixed-mode fracture criterion 

interfacial tractions 

nominal peel and shear strengths 

reduced peel and shear strengths 


Introduction 

Delamination is one of the predominant modes of failure in laminated composites when there is no 
reinforcement in the thickness direction. In order to design structures that are flaw- and damage-tolerant, 
analyses must be conducted to study the propagation of delaminations. The calculation of delaminations 
can be performed using cohesive elements 1 ' 2 , which combine aspects of strength-based analysis to predict 
the onset of damage at the interface, and fracture mechanics to predict the propagation of a delamination. 

Cohesive elements are normally designed to represent the separation at the zero-thickness interface 
between layers of three-dimensional elements (3D). Unfortunately, the softening laws that govern the 
debonding of the plies often result in extremely slow convergence rates for the solution procedure, and 
require very refined meshes in the region ahead of the crack tip. In addition, it is often difficult to model 
zero-thickness interfaces using conventional finite element pre-processors. 

The typical requirement of extremely fine meshes ahead of the crack tip 1 ' 3 ' 4 has been recently relaxed 
by a procedure that artificially reduces the interfacial strength while keeping the fracture energy constant 
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in order to increase the length of the cohesive zone. A longer cohesive zone allows the use of larger 
cohesive elements 3 . However, the difficulties associated with the generation of complex meshes 
composed of stacks of 3D elements joined together with zero-thickness cohesive elements often render 
the use of cohesive elements impractical outside research environments. 

Shell elements are more efficient for modeling thin-walled structures than 3D elements. Shell finite 
element analysis has been used in conjunction with the Virtual Crack Closure Technique (VCCT) by 
Wang et al. 5, 6 to evaluate strain energy release rates for the damage tolerance analysis of skin -stiffener 
interfaces. Wang used a wall offset to move the nodes from the reference surfaces to a coincident location 
on the interface between the skin and the flange. Therefore, compatibility was ensured by equating the 
three translational degrees of freedom of both shells with constraint equations and by allowing the 
rotations to be free. These models were shown to require a fraction of the degrees of freedom than are 
needed for full 3D analyses while producing accurate values for mode I and mode II strain energy release 
rates. 

However, Glaessgen et al. 7 observed that shell VCCT models of debond specimens in which the 
adherends are of different thickness predict the correct total energy release rate but that the mode mixity 
does not converge with mesh refinement. This phenomenon is similar to the well known oscillatory 
behavior of the in-plane linear elastic stress field near the tip of a bi-material interfacial crack 8 . 

Very few references discuss the use of shell models in combination with cohesive elements and most of 
them are designed for explicit time integration. A shell model for discrete delaminations under mode I 
loading using finite-thickness volumetric elements connecting shell elements was proposed by Reedy 9 . 
Johnson et al. 10 applied the bilinear traction-displacement of the Crisfield 1 1 cohesive model in conjunction 
with the Ladeveze mesomodel 12 for intralaminar damage to investigate wing leading edge impact and 
crash response of composite structures. A mixed-mode adhesive contact shell element that accounts for 
the thickness offset and the rotational degrees of freedom in the shells was developed by Borg 13 for the 
explicit finite element code LS-DYNA". Bruno et al. 14 developed shell models that use rigid li nks to 
offset the location of the cohesive element nodes from the shell reference surface to the interface between 
the layers. Using this technique, multi-layered shell models were analyzed using the zero-thickness 
cohesive elements available in the commercial finite element code LUSAS®. 

Cohesive elements have been found useful to study fracture along bi-metallic interfaces 15 " 18 . Cohesive 
laws eliminate the stress singularity at the delamination front and, consequently, they prevent the 
oscillation in the stress field and the corresponding non-convergence of fracture mode mixity 18 . 
Therefore, it is expected that the non-convergence of the mode mixity observed by Glaessgen in shell 
VCCT models, which is also caused by a mismatch of the elastic stiffnesses of the adherends, is also 
avoided in cohesive models. However, a study of the mesh independence of mode mixity in cohesive 
shell models has not been attempted by the present authors. 

The objective of this work is to present a method for the construction of shell cohesive elements for 
implicit analysis. The formulation of an 8-node cohesive element for modeling delamination between 
shell elements on the basis of a 3D cohesive element previously developed by the authors 1, 2 is proposed. 
The formulation takes into account the rotational degrees of freedom at the shell nodes and the relative 
location of the interface which allows the construction of models with more than one layer of cohesive 
interfaces and does not require coincident nodes along an interface nor does it require the use of rigid 
links. Therefore, the formulation allows the use of several layers of shell/cohesive elements through the 
thickness of a thick laminate. Since the cohesive elements have the topology of a standard 8-node 
element, the models are simpler to generate than those based on zero-thickness interfaces. 

In the following sections, the constitutive damage model previously developed by the authors will be 
outlined. Then, a procedure to construct a shell cohesive element from that of any conventional zero- 
thickness cohesive element will be described. Finally, three different examples of increasing structural 
complexity will be presented to show that simple shell models based on the proposed cohesive element 
for shells can be used to represent the onset and propagation of delamination in structural components. 
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Bilinear Cohesive Law 

A cohesive constitutive law relates the traction, x, to the displacement jumps, A, at the interface. The 
bilinear softening model, which is chosen here for its simplicity, is shown in Fig. 1 . One characteristic of 
all softening models is that the cohesive zone can still transfer load after the onset of damage (A n in Fig. 
1). For pure mode I, II or III loading, after the interfacial normal or shear tractions attain their respective 
interlaminar tensile or shear strengths, the stiffnesses are gradually reduced to zero. The area under the 
traction-displacement jump curves is the respective (mode I, II or III) fracture energy. Using the 
definition of the J integral, it can be shown that for small cohesive zones 19 : 

| q r (A) dA. = G c [1] 

where G c , the total area under the traction-displacement law, is the critical energy release rate for a 
particular mode, and A 1 is the displacement jump at failure, shown in Fig. 1. The penalty stiffness K is an 
arbitrarily large number selected such that the presence of undamaged cohesive elements does not 
introduce appreciable compliance to the structure 3 . 

A bilinear cohesive law is assumed for each of the three fracture modes. It is assumed that the 3 
direction is normal to the interface and that the interlaminar shear strength x° hear is independent of the 
shearing direction. Then, the displacement jumps for damage initiation in each mode are simply 
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similarly, the final displacement jumps are proportional to their corresponding toughness G ic as 
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Figure 1. 



4 


Cohesive Law for Mixed-Mode Delamination 

In structural applications of composites, delamination growth is likely to occur under mixed-mode 
loading. Therefore, a general formulation for cohesive elements must deal with mixed-mode delamination 
growth problems. 

Under pure mode I, II or III loading, the onset of damage at the interface can be determined simply by 
comparing the tractions with their respective allowable. Under mixed-mode loading, however, damage 
onset may occur before any of the traction components involved reach their respective allowable. 
Therefore, a mixed-mode criterion must be established in terms of an interaction between components of 
the energy release rate. 

The constitutive damage model used here, formulated in the context of Damage Mechanics (DM), was 
previously proposed by the authors 2 . The damage model simulates delamination onset and delamination 
propagation using a single scalar variable, d, to track the damage at the interface under general loading 
conditions. An initiation criterion that evolves from the Benzeggagh-Kenane fracture criterion 20 (B-K) 
ensures that the model accounts for changes in the loading mode in a thermodynamically consistent way 
and avoids restoration of the cohesive state. 

The constitutive model prevents interpenetration of the faces of the crack during closing, and a Fracture 
Mechanics-based criterion is used to predict crack propagation. The parameter A is the norm of the 
displacement jump tensor (also called equivalent displacement jump norm), and it is used to compare 
different stages of the displacement jump state so that it is possible to define such concepts as ‘loading’, 
‘unloading’ and ‘reloading’. The equivalent displacement jump is a non-negative and continuous 
function, defined as: 


A = J(A 3 ) 2 + (A s/war ) 2 [4] 

where <•> is the MacAuley bracket defined as <x>=4(x + |x|), which sets any negative values to zero. 

The term A 3 is the displacement jump in mode I, i.e., normal to the mid-plane, and A sW is the tangential 
displacement calculated as the Euclidean norm of the displacement jump in mode II and in mode III: 

a ^=a/(A i) 2 +(A 2 ) 2 [5] 

A bilinear cohesive law for mixed-mode delamination can be constructed by determining the initial 
damage threshold A° from the criterion for damage initiation and the final displacement jump. A 1 , from 
the formulation of the propagation surface or propagation criterion 2 . For the B-K fracture criterion, the 
mixed-mode displacement jump for damage initiation is 2 


A° = 




G„ + G w 


[ 6 ] 


where the B-K parameter r/ is obtained by curve-fitting the fracture toughness of mixed-mode tests and 
the fracture mode ratio is 2 


' G u + G u A P 2 

v g t ) \ + ip 2 -ip 


[ 7 ] 


where the displacement jump ratio is defined as P = A shear /(A shear + (A 3 )). 
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The displacement jump for final fracture is also obtained from the critical displacement jumps as 


a° 3 a{ +(a° w a{ w -a° 3 a{): 


r G ]l +G nL y 


tf=- 


[ 8 ] 


During overload, the state of damage d is a function of the current equivalent displacement jump A : 


A f (A-A°) 
(a" -A 0 ) 


The corresponding tractions can be written as 
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[ 10 ] 


where the Kronecker 8y is used to prevent the interpenetration of the surfaces of a damaged element 
when contact occurs 2 . 

In summary, by assuming that toughnesses and strengths in modes II and III are equal to each other, it 
is shown that the mixed-mode constitutive equations of a cohesive element are defined by five material 
properties and three displacement jumps: G Ic , G nc , P , r ( ’ hear , K, A,, A 2 , and A, . This formulation of the 
damage model allows an explicit integration of the constitutive model and ensures consistency in the 
evolution of damage during loading, unloading, and changes in mode mixity. The derivation and 
implementation of the cohesive damage model is described in Refs. [2, 21] and the reader is referred to 
those sources for details that will not be repeated here. 


Viscoelastic regularization 

Material models exhibiting softening behavior and shaip stiffness degradation typically exhibit severe 
convergence difficulties in implicit analysis programs. To improve the convergence rate of the iterative 
procedure, a viscous regularization scheme was implemented as suggested by ABAQUS 22 for its cohesive 
element and composite damage models. The procedure is a generalization of the Duvaut-Lions 
regularization model which, for sufficiently small time increments, maintains the tangent stiffness tensor 
of the softening material as positive definite. Defining d as the “kinematic” damage that depends on the 
displacement variables and d v as the “regularized” damage variable, the rate of change in the damage 
variable d v is a function of the damping factor p: 


d v = — (d — d v ) 

P 

The regularized damage variable at time t+At is updated as 

d v 


- & d 

1 P H 1 ’ 

t+At P +M l+At 

+ P+ti ° , 


[ 11 ] 


[ 12 ] 


It can be observed that the viscoelastic factor has the units of time and, therefore, its effect depends on 
the pseudo-time load proportionality factor which is selected arbitrarily by the user. In the following 
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sections, the effect of the pseudo-time on the viscoelastic factor is eliminated by specifying it as a 
nondimensional value, plt F , where t F is the estimated pseudo-time at the maximum load. 

The response of the damaged material is evaluated using the regularized damage variables. Using 
viscous regularization with a viscosity parameter that is much smaller than the time increment can 
improve the rate of convergence in the presence of softening processes without compromising the 
accuracy of the results. However, the use of larger values of the regularization parameter can cause a 
toughening of the material response and produce a stable propagation of delamination in problems with 
unstable delamination growth. 

Finite Element Implementation of Zero-Thickness 3D Element 


Zero-thickness cohesive elements as shown in Fig. 
connecting two laminae of a composite laminate. 


Cohesive element 
between solids 


Solid Element 



3 DOF per node 

Figure 3. a) 3D cohesive elements 


3 a are used to simulate the response of the interface 
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Shell Element 
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6 DOF per node 

b) shell cohesive elements. 


The constitutive equation of zero-thickness cohesive elements is established in terms of displacement 
jumps and tractions across the interface. The stiffness matrix for conventional 3D cohesive elements is 
constructed using standard isoparametric linear Lagrangian interpolation functions. In a local coordinate 
frame with x-y-z aligned such that the z-direction corresponds to the shell normal, the displacement jumps 
between the top and the bottom faces of the element are 


a 3 
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W 

A! 

> = • 

U 

, — 

U 

^2. 

3 D 
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top 
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= B U 


3D 


[13] 


bot 


where B is the matrix relating the element’s degrees of freedom U 3/ > to the displacements between the top 
and bottom interfaces. The Principle of Virtual Work for 3D cohesive elements can be written as 23 


S\}\ D B r ((i - D)c)B dA U 3D = SV T iD j ( B r S dA 


[14] 
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where D is a diagonal matrix composed the damage state d at the interface. The term I represents the 
identity matrix, and C is the undamaged constitutive matrix whose diagonal is composed of the penalty 
stiffnesses K, and S is the stress tensor. The stiffness matrix K 3D of the element is easily identified in the 
left hand side of Eq. 14 as the integral over the area A of the element: 

K 3 WX((l-D)c)BdT t 15 ] 

and the right hand side of Eq. 14 provides the element’s internal force: 

*, D =\ySdA [16] 

Because the stiffness K 3/J is constructed from the displacements U 30 , only three degrees of freedom can 
be active at a node, i.e., u, v, and w. The integration of the terms in Eqs. 15 and 16 is performed 
numerically using either Newton-Cotes or Gaussian integration. The element proposed is implemented in 
ABAQUS 22 as a user-defined element. 

Formulation Shell Cohesive Element 

Shell cohesive elements have a finite non-zero thickness because they connect nodes on the reference 
surfaces of two shells, as shown in Fig. 3b. In addition to the u, v, and w translational degrees of freedom, 
the nodes of a shell element also possess three nodal rotations 9 X , O y , and 0 Z . For kinematic continuity 
before delamination propagation, all three translational degrees of freedom at the interface must be the 
same in both adjacent plies, as shown in Fig. 4. 



layer 

interface 


Figure 4, Kinematic continuity is achieved when the displacements at the interface are equal in 
both plies 24 . 


When using shell elements based on First Order Shear Deformation, the three displacements at the 


interface are calculated from the nodal degrees of freedom as 23 : 
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Shell 


where tj is the thickness of the shell at node j, and its sign depends on whether the shell node described 
the displacements of the upper or the lower ply. For arbitrary orientations of a shell in space, the 
kinematic transformation matrix R k| is rotated such that the z-direction is aligned with the normal 



direction defined by two connected shell nodes. A transformation matrix R for the entire set of degrees of 
freedom in an element can be assembled using the transformations Rkj such that 

U 3 D=RU 5/ie „ and SV 3D =RSV shell [18] 

The stiffness matrix K mi of a cohesive element for shell analysis can be constructed by substituting 
Eqs. 18 into the expression of the Principle of Virtual Work for a 3D cohesive element (Eq. 12): 

<?uL//R r JX((i-d)c)b«14 ru„ =<?uL//R r IX SdA [19] 

By comparing Eqs. 14 and 19, it is clear that the stiffness and internal force of a shell cohesive element 
can be constructed from any standard cohesive element formulation using the following expressions: 

R Shell = R R 30 R an( l ^ Shell = R U III [20] 

Eqs. 20 show that a cohesive element for shells is easily constructed from any 3D cohesive element by 
simple matrix multiplication. An 8-node cohesive element for shells was implemented as an ABAQUS 
UEL (user element) subroutine 22 by using a zero-thickness 3D element previously developed by the 
authors 2 as a kernel. 


Spatial Orientation of the Shell Normals 

The transformations shown in Eqs. 20 is only valid when the normals to the upper and the lower shells 
are aligned with the z-direction of the global reference frame. In the case of arbitrary orientation of the 
shells in space, it is necessary to define local reference frames at each comer of the shell cohesive 
elements. In the present model, the local normal direction (z') is defined by the vector between the lower 
and upper nodes at a comer. Using a rotation tensor, the normal and tangential components of the 
displacement jump can be expressed in terms of the displacement field in global coordinates. 


Special Configurations 

It was shown in Eqs. 17 and 20 that the element formulation is based on the thickness t" and t L of the 
upper and lower shells being connected together by a cohesive element. Therefore, the values of the 
thicknesses must be provided for each comer of every cohesive element in the model by using element 
property data cards. In models where the thicknesses of the adherends vary across the model, different 
property cards would have to be generated for every cohesive element. To improve the ease of generating 
a model where the shell thicknesses are non-uniform, it is convenient to define common geometric 
configurations that allow the use of single property cards for large areas of cohesive elements. Two such 
configurations are shown in Fig. 5. The configuration shown in Fig. 5a, which is typical of skin/flange 
constructions, can be defined by specifying the thickness f of the skin, which is constant for all cohesive 
elements in the feature. For configurations that are symmetric with respect to the interface, as in the 
example shown in Fig. 5b, the ratio t u /t l is constant. In either case, the two thickness t " and t l can be 
determined within the element with the additional knowledge that the distance t“ + t l is equal to twice the 
distance between the upper and the lower nodes and a single property card can be used for all elements in 
the feature. 
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a) t l = constant 



Figure 5. Special configurations with constant thickness properties. 


Example and Verification Problems 


Problem 1: Mixed Mode Bending (MMB) Specimen 

The MMB test offers the possibility of measuring different mixed-mode ratios, ranging from pure 
Modes I to II, by changing the length e of the loading lever shown in Fig. 6 . The loading arm is modeled 
with a multi-point constraint equation between the displacements at the load application point and the two 
loading points on the upper arm of the specimen. 

The MMB test of a 24-ply unidirectional AS4/PEEK (APC2) carbon fiber reinforced composite was 
simulated with a simple model composed of two shell layers, and the results were compared with the 
results of a 3D model previously published by the authors 2 . The specimen is 102 mm-long, 25.4mm- 
wide, with two 1.56 mm-thick arms. The initial delamination length is a 0 =34mm, and the length e of the 
lever is 43 mm, which was calculated to provide a mode mixity G/Gf= 50%. The material properties of 
the specimen are shown in Table 1, and the properties of the interface are shown in Table 2. It should be 
noted that the parameter 77 was calculated from a fit of experimental data corresponding to a single test 
specimen at each mode ratio 25 . 



Figure 6. MMB test configuration and finite element models. 
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Table 1. Material Properties of AS4/Peek ?5 . 


En(GPa) E 22 -E 3 3 (GPa) v 12 =vi 3 V 23 Gi 2 -Gi 3 (GPa) G 2 3 (GPa) 

122.7 10.1 0.25 0.45 5.5 3.7 


Table 2. Properties of the interface 1 ' 25 . 

G Ic (N/mm) G IIc (N/mm) r 3 ° (MPa) r° w (MPa) // 

0.969 1.719 80 100 2.28 


The experimental results shown in Fig. 7 were performed by Reeder 25 . They relate the load to the 
displacement of the point of application of the load F on the lever. The analytical solutions shown in the 
figure were calculated using the equations from Reeder 26 with G c =1.123 N/mm, which corresponds to the 
point where the delamination growth was first observed. This point is the second experimental data point 
in Fig. 7. Therefore, this point lies on the analytical constant-G curve. It can be observed that the test data 
indicates an R-curve effect consisting of an increase in toughness that ranges from an initiation value of 
1.123 N/mm to a propagation value of approximately 1.27 N/mm. 

The shell finite element model shown in Fig. 6 was developed using two layers of shell elements. Each 
layer is composed of a mesh of 2X 150 elements. Therefore, the size of the elements in the direction of 
crack propagation is 0.68mm. The shell model only represents 1/10 of the specimen width, so the 
predicted applied load was multiplied by a factor of 10. The loading lever was modeled as a rigid link 
using a multi-point equation. To improve the convergence rate, a viscoelastic factor of plt F = 2.0-10~ 4 

was used, where t F is the estimated pseudo-time at the maximum load. A penalty stiffness of K= 10 5 
N/mm 3 was used. Recommendations on the determination of K are provided in Ref. 4. The corresponding 
load-displacement results are shown in Fig. 7. For comparison, the results of a simulation using a three- 
dimensional model by the authors 2 is also shown in Fig. 7. It can be observed that both models predict the 
analytical peak load and load-displacement curve accurately. The correlation of the analytical and finite 
element models with the test results could be improved using the toughness for propagation rather than 
the value for initiation. For an applied displacement larger than 7 mm, the delamination front proceeds 
under the mid-span load point and the present analytical solution is no longer valid. At that point, the 
finite element models diverge from the analytical solution. 

To verify the effects of mesh size and viscoelastic regularization, a new shell model was created with 
300 elements in the propagation direction, or twice the number of elements of the base model. The model 
was analyzed with p=0. The predicted load with the refined mesh model is less than 1% lower than the 
load obtained from the previously discussed coarser model with pi t F = 2. • 1(T 4 , as can be observed in Fig. 
8. Flowever, the refined mesh model failed to converge after a short amount of crack propagation. 

An alternative to using fine meshes of cohesive elements was proposed by Turon et al. 3 The approach 
consists of reducing the strengths and r° shear such that the length of the cohesive zone increases to 
cover at least three elements. According to Turon, the strengths that are required to achieve convergence 
using an element size of 0.68 mm are: 


* 



J_ l 9*£ 33 G Ic _ Q 
r 3 ° ]j 32x5/ 


^ shear _ ^ 1 9/T E^ 3 Gjj c _ ^ 

Ahear Ahear V 32x5I eIe 


[ 21 ] 


Equations 21 indicate that the strength reduction needed for mesh convergence is greater for peel than 
for shear. Since both strengths should be reduced by the same factor to maintain a constant mode mixity, 
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the reduced strengths r* and T* hear that ensure that the cohesive zone is covered by at least five elements 
are: = 0.63 xr 3 ° = 50.4 MPa ; T* hear = 0.63 xr s 0 tear = 63.0 MPa . 



Applied displacement, mm. 

Figure 7. Predicted and experimental response of 50% mode mix MMB test. 



Figure 8. Effect of mesh density and viscoelastic regularization on predicted load displacement 
response of 50% mode mix MMB test. 


12 


An analysis was performed with the original (coarse) mesh model with reduced strengths and no 
viscoelastic regularization. The results shown in Fig. 8 indicate that the results from the original (coarser) 
mesh using reduced strengths are virtually identical to those using a refined mesh. In the present problem, 
the use of lowered strengths reduces the occurrence of divergence in the iterative solution and the use of 
viscoelastic regularization is not required. 


Problem 2: Debonding of Skin-Stiffener Specimen 

Most composite components in aerospace structures are made of panels with co-cured or adhesively 
bonded frames and stiffeners. Testing of stiffened panels has shown that bond failure at the tip of the 
stiffener flange is an important and very likely failure mode. Comparatively simple specimens consisting 
of a stringer flange bonded onto a skin have been developed by Minguet et al. 27 to study skin/stiffener 
debonding. The configuration of the specimens shown in Fig. 9 was analyzed by Cvitkovich et al. 28 The 
specimens are 203 mm-long, 25.4 mm-wide. Both skin and flange were made from IM6/3501-6 
graphite/epoxy prepreg tape with a nominal ply thickness of 0.188 mm. The skin lay-up consisting of 14 
plies was [0/45/90/-45/45/-45/0]s, and the flange lay-up consisting of 10 plies was [45/90/-45/0/90]s. 

The properties of the unidirectional IM6/3501-6 unidirectional graphite/epoxy tape reported by 
Cvitkovich 28 are shown in Table 3. However, the interface properties for IM6/3501-6 were not available. 
Instead, the interface properties of AS4/3501-6 are used, as suggested by Cvitkovich 28 . The values of Gi c , 
Gji c and the parameter 7=1.75 for the B-K criterion shown in Table 4 are taken from test data for 
AS4/3501-6 29 . 


Tension Test 


203 mm 


Strain gages 



42 mm 
50 mm 


127 mm 

Extensometer 

Figure 9. Skin-stiffener specimen configuration 28 . 

Table 3. Material Properties of IM6/3501-6 unidirectional graphite/epoxy tape [Ref. 28] 


E n (GPa) E 22 -E 33 (GPa) v 12 =vi 3 

v 23 

G 12 =G 13 (GPa) 

G 23 (GPa) 

144.7 9.65 0.3 

0.45 

5.2 

3.4 


Table 4. Properties of the interface [Refs. 29, 28, 30]. 

G lc (N/mm) G IIc (N/mm) f (MPa) f hear (MPa) 7 
0.0816 0.554 61 68 1.75 
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A complete analysis of the delamination growth in the tension specimen requires a high degree of 
complexity. For instance, Krueger 31 developed highly detailed two-dimensional models using up to four 
elements per ply thickness. In previously published work 30 , the authors determined that it is possible to 
predict the debond load of the specimen using cohesive elements using a relatively coarse three- 
dimensional model. To keep the modeling difficulties low and the approach applicable to larger problems, 
the 3D model was developed using only two brick elements through the thickness of the skin, and another 
two through the flange. The complete model consists of 1,002 three-dimensional 8-node C3DI elements 
and 15,212 degrees of freedom. This model did not contain any pre-existing delaminations. 

To account for the effect of residual thermal strains, a thermal analysis step was conducted before the 
mechanical loads were applied. The same coefficients of thermal expansion («//=-2.4x 1 0-8/°C and 
« 22 =3.7x 1 0-5/°C) are applied to the skin and the flange, and the temperature difference between the 
curing and room temperature is AT= -157 °C [Krueger 31 ]. The flange has more 90° plies than 0° plies and 
the skin is quasi-orthotropic, so that residual thermal stresses develop at room temperature at the 
skin/flange interface. 

In contrast to previous models of the skin -stiffener specimen, an extremely simple shell model was 
developed that uses one layer of shell elements to represent the skin and another layer to represent the 
flange of the stiffener. When subjected to tensile loads, the primary mechanism for debonding is shear 
lag. This fracture mechanism is captured in the model by using a layer of shell cohesive elements to bond 
the two shell layers together. The model is composed of six elements across the width, and the mesh in 
one of the two flange tip regions is refined. The element size in the propagation direction is 0.5 mm, and 
the rest of the flange is meshed with 1.0 mm elements. A penalty stiffness of K= 10 5 N/mm 3 was used for 
this analysis. The model is composed of 624 shell and cohesive elements and only 3192 degrees of 
freedom. Since this particular problem is dominated by the initiation of the debond rather than its 
propagation, the nominal strengths were not reduced, and an accurate calculation of the energy release 
rates during propagation was not attempted. 

A deformed plot of the finite element model immediately after flange separation is shown in Fig. 10. 



Figure 10. Two-shell model of skin-stiffener specimen. 

In previous work by the authors, it was stated that eigenmode analysis of the element stiffness matrices 
has shown that Gaussian integration can cause undesired spurious oscillations of the traction field when 
large traction gradients are present over an element and, consequently, that Newton-Cotes integration is 
superior to Gaussian integration in cohesive elements 1 ' 32 . However, in the analyses of the skin-stringer 
specimen, it was found that the Newton-Cotes integration causes odd wrinkling of the thin flange tip 
during the onset of delamination, as illustrated in Fig. 11a. The wrinkling deformations were found to 
delay the full separation of the flange tip. On the other hand, the wrinkling deformations were virtually 
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eliminated by using Gaussian integration, as shown in Fig. lib. Consequently, previously established 
conclusions on the preferred integration scheme may need to be re-examined. 




a) Newton-Cotes integration 


b) Gaussian integration 


Figure 11. Onset of delamination with two different integration schemes (10X deformation). 


The predicted load-extensometer response for two values of viscoelastic regularization are shown in 
Fig. 12. The initiation of delamination predicted using a viscoelastic regularization factor of 
plt F = 2. • 10 4 is 17.2 kN, which is 24% below the mean debond load measured by Cvitkovich 28 (22.7 
kN). The unstable nature of delamination propagation can be deduced from the reduction of the applied 
load during the propagation of the delamination. The predicted delamination initiation load using 
ptt F = 2.-KT 3 is 17.5 kN. However, delamination propagation is stabilized by the use of a larger amount 
of regularization, as can be observed in Fig. 12. The advantage of using a higher amount of regularization 
is that the solution is obtained much faster than with a lesser amount of regularization: the model with a 
regularization factor of 2. TO -4 took 6 minutes of CPU time on a 3.20GHz Xeon-based PC, while the 
model with a factor of 2.- 10 3 took less than 3 minutes. 

The previously published results of a 3D model are also shown in Fig. 12. The slope of the force- 
extensometer curves for the 3D and shell models are different because nodal displacements were used to 
calculate the extensometer measurement. For the 3D model, the nodes are located on the surface of the 
laminates, while for the shell models, the nodes are located on the reference surfaces. It can be observed 
that the 3D model predicts a more accurate debond load of 21.9 kN, which is probably a result of the 
higher fidelity of the flange tip region in the 3D model. 
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Figure 12. Predicted and experimental load-extensometer response for skin-stiffener specimen in 
tension. 


Problem 3: Thick Composite Lug 

Shell finite elements are typically reserved for the analysis of thin-section structures such as fuselage 
skins. Since shell elements either neglect or simplify the transverse deformations, they are not usually 
applied to model components in which the thickness dimension is of the order of magnitude of other 
characteristic dimensions of the component. However, computationally efficient models of thick 
components have occasionally been developed by stacking layers of shell elements. For instance, the 
complex through-the-thickness deformations of thick sections of the Composite Armored Vehicle, which 
includes layers of glass/epoxy, ceramic tiles, and rubber, were successfully modeled by tying four or more 
layers of shell elements using multi-point constraint equations 24 . 

In addition to the computational efficiency that results from using shell elements, the use of shell 
cohesive elements to model layered structures provides the added benefit of simplifying the construction 
of complex models in several ways. Firstly, shell cohesive elements eliminate the laborious task of 
constructing multi-point constraints to tie the layers together. Secondly, unlike conventional cohesive 
elements, shell cohesive elements have a non-zero thickness and they have the same topology and 
connectivity as a standard volume element. Therefore, they can be constructed using any finite element 
pre-processor. Finally, shell cohesion elements do not require the use of wall offsets. Therefore, the 
construction of shell models is simplified and, more importantly, the shell reference surfaces do not have 
to be shifted to a co mm on interface, which allows stacks of multiple shell layers. 

During the course of a failure analysis of a composite aircraft fin undertaken at NASA Langley 33 , a 
layered shell element of the critical attachment lug area was developed. The model detail shown in Fig. 
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13 illustrates how cohesive elements for shells can be used to construct models of large and geometrically 
complex thick composite structures by stacking layers of shell elements. The thick lug fittings were 
modeled with 14 layers of shell elements, which were connected with shell cohesive elements. All other 
regions of the model were modeled with a single layer of shell elements. The pin assembly was modeled 
as a rigid surface with a diameter equal to that of the lug hole. Frictionless contact equations were 
prescribed between the edge of the shells around the bolt hole and the rigid surface. This model has 20886 
nodes with 34524 elements. The response of the layered-shell model was compared to a 3D model (no 
delamination) and it was found that the strain distributions predicted by the two models were in close 
agreement 31 . 



Figure 13. Layered shell model of composite lug with ply damage and cohesive elements. 


In addition to delamination at any interface between the 14 layers of shells, the model uses a 
progressive damage model in which the matrix and fiber damage is simulated by degrading the material 
properties. An example of damage analysis is shown in Fig. 13, where the regions in red represent areas 
where damage is predicted. The analyses for intra-ply damage were developed within ABAQUS using a 
USDFLD user- written subroutine. The finite element program calls this subroutine at all integration 
points of elements that have material properties defined in terms of the field variables. The subroutine call 
provides access points to a number of variables such as stresses, strains, material orientation, current load 
step, and material name, all of which can be used to compute the field variables. Stresses and strains are 
calculated at each incremental load step, and evaluated by the failure criteria to determine the occurrence 
of failure and the mode of failure. 

The results of a failure analysis indicate that a cleavage-type failure of the lug, in which an initial 
fracture develops in a direction parallel to the orientation with the greatest number of fibers. This 
direction is nearly aligned with the load direction, as shown in Figs. 14. 
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Figure 14. Progressive failure analysis predicts a two-step cleavage failure of the composite lug 33 . 

The greatest advantage of the layered shell model is its ability to predict the initiation of delamination, 
as well as the component’s tolerance to delamination flaws. A study was undertaken to evaluate the 
damage tolerance of the lug in the presence of a pre-existing delamination. The flaw was introduced in the 
regions of the model shown in red in Figs. 15 by setting the region of shell cohesive elements shown in 
red to “failed” by defining the damage variable of these elements as unity. Failed cohesive elements do 
not have stiffness in tension or shear, but they prevent interpenetration of the crack surfaces in 
compression. 



Figure 15. Configuration of delaminated lug and model for damage tolerance analysis. 

The strength analyses with and without a pre-existing delamination were performed under tensile 
loading. For this loading condition, the analyses did not predict delamination initiation nor propagation of 
the delaminated region. The load-displacement results of the two analyses are shown in Fig. 16. It can be 
observed that the initial flaw does not appreciably affect the overall stiffness or strength of the composite 
lug. The predicted loads for both simulations fall between the load measured for the subcomponent test 
performed in 2003 and the W375 tension load that is believed to have caused the in-flight failure of the 
lug 33 . 
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1200 


1000 


Subcomponent test, 2003 



Figure 16. Predicted lug strength with and without initial flaw. 


Concluding Remarks 

The formulation of a cohesive element for shell analysis that can be used to simulate the initiation and 
growth of delaminations between stacked, non-coincident layers of shell elements was presented. The 
procedure to construct the element accounts for the thickness offset by applying the kinematic relations of 
shell deformation to transform the stiffness and internal force of a zero-thickness cohesive element such 
that interfacial continuity between shell layers is enforced. The internal load vector and the stiffness 
matrix of a cohesive element connecting shells can be calculated in a very simple way by means of a 
transformation matrix, independently of the cohesive constitutive models. 

Stacks of shell elements can be used to create effective models that predict the inplane and 
delamination failure modes of thick components. The formulation allows the use of multiple layers of 
shells through the thickness and it does not require offsetting nodes to a common reference surface. In 
addition, shell cohesive elements have the same topology as three-dimensional elements so they can be 
generated using any finite element pre-processor while zero-thickness cohesive elements for three- 
dimensional elements cannot be generated with most pre-processors. 

Using representative examples of increasing structural complexity, it was shown that the proposed 
cohesive element for shells can be used to represent the propagation of a pre-existing delamination, as 
well as the onset and propagation of delamination in components that do not contain pre-existing cracks. 
The analysis of a 50% mode-mix MMB test indicates that the initiation and propagation of delamination 
under mixed mode can be predicted accurately. It was found that the convergence difficulties that are 
typical of cohesive laws can be mitigated with the use of viscoelastic regularization or by reducing the 
interfacial strength to enlarge the process zone. The analysis results of a two-layer model of a skin- 
stiffener debond indicate that such a simplified model can produce results that may be sufficiently 
accurate for design and down-select analyses. For better predictive accuracy, models in which the flange 
taper is represented more accurately would be required. Finally, the analysis results of a large fin lug 
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component was presented to illustrate how layered shell models with cohesive elements can be used to 
predict the strength of thick composite structures. This combination of cohesive elements and shells 
represents an effective way of simulating complex aeronautical structures that incorporate other non- 
linearities such as ply damage propagation, contact, and finite rotations. However, for the loading 
conditions that were analyzed, the cohesive elements in the model were not subjected to loads sufficiently 
large to cause any delamination growth. 
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